// program to calculate skewness from delta files; one can generate delta files from IC's as well if one wants to find skewness directly from ICs.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>


FILE *fopen_check(char *filename, char *mode)
{
  FILE *fp = fopen(filename, mode);
  if (fp == NULL) {
    fprintf(stderr,"Error opening file %s with mode %s\n",filename,mode);
    exit(1);
  }
  return fp;
}

int main(int argc, char *argv[]) {
  int ng=128, ng3, i;
  FILE *fp;
  
  fp = fopen_check(argv[1], "rb");
  fseek(fp, 0, SEEK_END);
  ng3=ftell(fp)/4;
  ng=rint(cbrt(ng3));
  if (ng*ng*ng != ng3) {
    fprintf(stderr, "Error: appear to be %d floats in %s, not a cube?\n",
	    ng3,argv[1]); exit(EXIT_FAILURE);
  }  
  
  float delall[ng3];
  rewind(fp);
  fread(delall, 4, ng3, fp);
  // now loop over all deltas and calculate the <d^3> and <d^2> to get M3
  double d1=0., d2=0., d3=0., d4=0., d5=0.;
  for (i=0; i<ng3; i++) {
    d1+=delall[i];
    d2+=delall[i]*delall[i];
    d3+=pow(delall[i], 3.0);
    d4+=pow(delall[i], 4.0);
    d5+=pow(delall[i], 5.0);
  }
  d1/=ng3; d2/=ng3; d3/=ng3; d4/=ng3; d5/=ng3;
  printf("1 %g\n2 %g\n3 %g\n4 %g\n5 %g\n", d1, d2, d3/pow(d2,1.5), d4/pow(d2,2.0),d5/pow(d2,2.5));
  fclose(fp);
  
}


